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Variational wave functions used in the variational Monte Carlo (VMC) method are exten- 
sively improved to overcome the biases coming from the assumed variational form of the wave 
functions. We construct a highly generalized variational form by introducing a large number 
of variational parameters to the Gutzwiller-Jastrow factor as well as to the one-body part. 
Moreover, the projection operator to restore the symmetry of the wave function is introduced. 
These improvements enable to treat fluctuations with long-ranged as well as short-ranged cor- 
relations. A highly generalized wave function is implemented by the PfafRans introduced by 
Bouchaud et al, together with the stochastic reconfiguration method introduced by Sorella for 
the parameter optimization. Our framework offers much higher accuracy for strongly correlated 
electron systems than the conventional variational Monte Carlo methods. 
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1. Introduction 

Strongly correlated electron systems have brought 
many fundamental and challenging issues in condensed 
matter physics.^ They are characterized by a competition 
between the itinerancy of electrons favored by the kinetic 
energy and the localization caused by the Coulomb inter- 
action. When the latter contribution becomes predomi- 
nant, the material turns from metal into the Mott insu- 
lator at specific electron densities.^ This metal-insulator 
transition is called the Mott transition. Mott insulators 
and its related materials show fruitful properties such 
as the high-temperature superconductivity in copper ox- 
ides.^ Such phenomena are certainly beyond the frame- 
work of the standard band theory based on the one-body 
approximation. Many-body correlation effects play cru- 
cial roles in the strongly correlated systems. 

Theoretical routes to investigate these systems are 
severely restricted because of difficulties in treat- 
ing strong correlation effects. For this purpose, 
there exist several numerical methods, such as the 
exact diagonalization (ED), auxiliary- field quantum 
Monte Carlo (AFQMC),^""^ density matrix renormal- 
ization group (DMRG),* dynamical mean- field the- 
ory (DMFT),^'^° path- integral renormalization group 
(PIRG),"-!*^ Gaussian-basis Monte Carlo (GBMC) 
and variational Monte Carlo (VMC)^° methods. 

Among them, the VMC method is tractable in rel- 
atively large system sizes even at large amplitude of 
interactions and geometrical frustrations. However, the 
bias inherently and inevitably contained in the assumed 
variational form of the wave functions is a fundamental 
drawback in the VMC method. Therefore, construction 
of highly accurate wave functions is crucially important. 
In an interesting region where various phases compete, 
wave functions which do not sufficiently take into account 
quantum fluctuation effects often give even qualitatively 



wrong results. 

The VMC method^° offers the exact treatment of 
Jastrow-type wave functions^^ within the statistical ac- 
curacy. The Gutzwiller-Jastrow factors, which are 
operated to one-body wave functions, enable to take ac- 
count of many-body correlation effects and go beyond 
mean-field descriptions. However, in conventional treat- 
ments, the one-body parts are usually simply taken as 
the ground state of the mean-field Hamiltonian where 
the symmetry is explicitly broken by the mean field. The 
limitations and drawbacks of these conventional varia- 
tional wave functions are the following: 

(i) It is hard to describe different competing phases 
within a single variational form. 

(ii) They do not often satisfy inherent symmetry prop- 
erties because of the symmetry-broken one-body 
part. 

(iii) Although the one-body part crucially determines 
fundamental properties of the variational wave 
function, the one-body part remains primitive if 
quantum fiuctuation effects are not taken into ac- 
count. 

Recently, numerical techniques to optimize a huge 
number of variational parameters in the VMC framework 
are developed. ^^"^^ These developments have opened the 
possibility of overcoming a biased nature of the varia- 
tional approach and allow us to extend the potential of 
variational wave function study. The biases inherently 
and inevitably contained in the assumed variational form 
of the wave functions are aimed to be largely relaxed 
by a large number of variational parameters, which al- 
low us to treat fiuctuations with long-ranged as well as 
short-ranged correlations. One of the successful results 
in this approach is seen in the electron-state calculations 
for small molecules. ^^"^^ In these studies, a linear combi- 
nation of multi-configurational Slater determinants with 



2 J. Phys. Soc. Jpn. 

the Jastrow factor is chosen as a variational wave func- 
tion. All the parameters, such as linear coefficients, or- 
bitals in Slater determinants, and Jastrow parameters, 
are optimized by the recently developed energy mini- 
mization techniques. These efforts offer a reliable method 
to obtain quantitatively accurate wave functions in small 
molecules. However, this treatment can not be directly 
applied to the bulk electron systems, because we have to 
deal with hundreds of electrons and can not handle multi- 
configurational Jastrow wave functions within practi- 
cal computational costs. In the VMC studies on lattice 
models, a large number of variational parameters have 
been introduced to the Jastrow factor by Sorella.^^ This 
improvement has opened possibility to describe quan- 
tum phase transitions within a single variational wave 
function. ^'^ However, biases coming from the one-body 
part still remain because the conventional one-body part 
corresponds to the mean-field single Slater determinant 
with only a few variational parameters. This hypothesis 
strongly influences variational results even though the 
Gutzwiller-Jastrow factor with many parameters is in- 
troduced. In order to reduce the conventional biases, the 
one-body part must be reconstructed by a deliberately 
examined parameterization. 

In this paper, we reconstruct and improve the one- 
body part by introducing many variational parame- 
ters with well-thought-out and computationally tractable 
forms. We also introduce several symmetry projections in 
the ground state. We demonstrate the efficiency and ac- 
curacy of our variational framework in which many varia- 
tional parameters and the symmetry projection allow re- 
ducing the biases and bring us quantitatively more accu- 
rate wave functions than those in the literature. Our goal 
is to introduce a conceptually new scheme for strongly 
correlated electrons under large quantum fluctuations. 
This is crucially important in simulating regions near the 
quantum critical points and regions of competing orders 
with enhanced fluctuations. 

Single-band Hubbard model is suited for benchmark 
of many-body correlation effects. We improve variational 
wave functions in order to study the ground-state prop- 
erties of the Hubbard model on a square lattice defined 
by 

n = ^e{k)cl^Cka + U^ni-(nii, (1) 

k,<T i 

where 

4. = 7^ E 4e"=-' , = ^ c,.e-"=- (2) 

are the creation and annihilation operators. The number 
operator is njo- = cJ^Cjcr- The energy dispersion e{k) is 

given by 

e{k) = — 2t(cos kx + cos ky) — At' cos k^ cos ky, (3) 

where t {t') is the transfer integral between the Wan- 
nier orbitals of nearest-neighbor (next nearest-neighbor) 
sites. Model parameters are t'/t, U/t, and filling n = 
N/Ns, where N is the total number of electrons and Ng 
is the total number of sites. We take Ng = L x L sites 
with the boundary condition periodic in x direction and 
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antiperiodic in y direction (periodic-antiperiodic bound- 
ary condition). 

The organization of this paper is as follows. In §2, we 
introduce variational wave functions used in this study. 
The functional form with a large number of variational 
parameters and the symmetry projection are described. 
The optimization method to efficiently handle many vari- 
ational parameters is explained in §3. The accuracy of 
our variational wave functions is benchmarked by com- 
parisons with results obtained from unbiased methods in 
§4. Section 5 is devoted to summary and discussions. 

2. Variational Wave Functions 

In this section, we introduce variational wave functions 
used in this paper. We construct wave functions which 
bear the following properties: 

(i) Flexibility to describe several different phases by 
controlling many variational parameters of a unified 
variational form, 

(ii) Capability of treating many-body correlations be- 
yond mean- field states, 

(iii) Conservation of symmetry with quantum numbers 
expected in the ground state. 

Recent development in the VMC method allows us to 
deal with a large number of parameters. These nu- 
merical techniques are described in §3. We construct 
wave functions which enhance the capability of removing 
biases posed on the variational form. This is achieved at 
least partially by introducing enormous number of pa- 
rameters. 

2. 1 Functional form of variational wave functions 

The general functional form of wave functions in this 
paper is 

\^)=vm, (4) 

where |(/)) is a Hartrec-Fock-Bogoliubov type wave func- 
tion called "one-body part," £ is the quantum-number 
projector^^' controlling symmetries of wave function, 
and V is the Gutzwiller-Jastrow factor^"'^'^^ including 
many-body correlations. In order to improve variational 
wave functions within the sector classified by quantum 
numbers, we only employ V that preserves symmetries 
of C\4>). This means that £ and V are commutable 
{VC = CV). 

2.1.1 One-body part 

The one-body part usually corresponds to the mean- 
field Slater determinant with several variational pa- 
rameters. Though the Gutzwiller-Jastrow factor intro- 
duces many-body correlations, this variational hypoth- 
esis strongly restricts flexibility of wave functions. We 
reexamine the functional form of the one-body part and 
introduce as many as possible variational parameters in 
order to improve wave functions. 

First, we consider a Hartree-Fock-Bogoliubov type 
wave function with antiferromagnetic (AF) and super- 
conducting (SC) orders which have been introduced 
by Giamarchi and Lhuillier."^^ Here, we start from a 
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slightly different representation introduced by Himeda 
and Ogata. ^'^ The wave function diagonalizes the mean- 
field Hamiltonian 



feeAFBZ 



+ e\k)(bi^bk^+bi^bki) 



(5) 



where AFBZ denotes the folded AF Brillouin zone, 
Z\g(-,(fc) and Z\g(^(fc) are SC order parameters, and /i 
is the chemical potential. The energy dispersion of AF 
bands £"(fe) and e^{k) are defined as 

p»W=6W-y!WT^ (..AFBZ) (6) 

with a(fe) = - + Q))/2 and = + 

e{k + Q))/2, where Z\af is the AF order parameter and 
the vector Q corresponds to the AF order {Q = (tt, tt)). 

The operators {0^0-' ^lai iWkcr,bkcr}) are creation (an- 
nihilation) operators for the AF quasiparticles related 
to the electron operator c^^ {ckcr) through the following 
unitary transformation 



fe e AFBZj (7) 



with 



Uk{Vk) = 



m 



1/2 



(8) 



The wave function with N particles extracted from the 
eigenfunction of Hmf is written as 



'AF+SC/ 

with 

ip^ik) 

/(fe) : 



^ (^nfe)at^aL 

.feeAFBZ 



N/2 



(9) 



^ic(fe) 



(£«(fc) - + v(e«(fe)-/^)^ + [4ic(fe)F 

{e\k) - M) + ^(£f(fe)-/x)2 + [Zi|c(fe)]^ 



(10) 

In the limit Z\af 0, the operators al.^ and 6^^,^ 
are reduced to c^^ and c^^g^, respectively. The wave 
function |</>af+sc) (cq- (9)) becomes the conventional 
BCS wave function. On the other hand, in the limit 

zigc(fe)(^sc(^)) ^ 0' <P"(fe)(<p''(fe)) goes to zero if 
e"'{k){e^{k)) > ji and otherwise diverges. Thus, the AF 
quasiparticles are only filled below the chemical poten- 
tial and I^AF+sc) is reduced to the normal AF mean- field 
wave function. 



In VMC studies on lattice systems, several variational 
parameters are considered to improve the one-body part. 
For example, Z\af and A'^q [k) are adopted to describe 
the magnetism and superconductivity, respectively. The 
chemical potential ^ and band renormalization effects^'^ 
improve the accuracy. These variational parameters allow 
optimizations of Mfe, Vk, 'P°'{k), and i/^''(fe) to a certain but 
restricted extent in eqs. (8) and (10). 

Now, instead of taking Aaf, '^sci^)' ^^'^ A* as vari- 
ational parameters, we take parameters Uk, Vk, <^"(fe), 
f''{k) independently for each fe e AFBZ under the con- 
ditions: 

W'fc + = 1 1 ""-fc = "fe 1 "-fe = «fe > (11) 
^'^(-fe) = ^"(fc), /(-fe) = /(fc). (12) 

From eqs. (7) and (9), |<?!)af-i-sc) can be written with c^^: 



PAF-I-SC 



.feeAFBZ 

+ 



^ \[ul^'^{k)-vl<p\k))clc[ 

;AFBZ 

[-vlv'^{k)+ulv\k))ci 

+ (ip''{k) + ip\k)^ukVk 

^ ('^fe-l-QT'^-fei ~ '^feT^-fe-Qi 



k+Qt^-k-Ql 



N/2 



|0). 
(13) 



(14) 



Then, we transform Uk, Vk, <p"(fe), '^*(fe) to 
Uk = cos 9k, Vk = sin 6k 
A{k) = cos2 6lfe(p»(fe) - sin^ 6lfe/(fe) 
B{k) = - sin^ Ok^p^ik) + cos^ ek<p''ik) 

The coefficient of the third term in eq. (13) is rewritten 

as 

[if^ik) + ip\k))ukVk = (A{k) + S(fc)) tan^fc = C(fc). 

(15) 

As a result, variational parameters Uk, Wfc, ip°'ik), 'p^{k) 
are mapped to new parameters A{k), B{k), C{k). The 
parameter A{k) {B{k)) corresponds to singlet pairings in 
(out) the AFBZ. Finally, with the definitions (^(^^(fe) = 
^(fe), 1^(1) (fe + Q) = B{k), and ip'^^^k) = C{k) for fe e 
AFBZ, we obtain the wave function 



E ^^'Hfe)4ci,, 



.fcgBZ 



-1 N/2 



4 



k+Q^^-ki 



CfeT^-fe-Qi; 



feeAFBZ 



with the conditions 

ip(^\-k) = (^(i)(fc) , ip'^^\-k) = ip^^^\k). 



|0) 
(16) 

(17) 



The AF mean-field state is realized by using the sec- 
ond term. Dealing with ip'^^\k) and (^(^^(fc) directly 
as fe-dependent variational parameters allows us to ex- 
press various states such as paramagnetic metals, an- 
tiferromagnetically ordered states, and superconducting 
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states with any gap function within a single framework 



of 



r). Moreover, since the number of the variational 



parameters increases scaled by the system size, it allows 
taking account of fluctuation effects with short-ranged 
correlations. In this paper, we call Impair) a "generalized 
pairing hmction" and </7*^^)(fc), </?*^^^(fc) are called "pair 
orbitals." Introducing all the possible ordered vectors Q 
would further generalize |^^>pair)- However, this extension 
substantially increases the number of variational param- 
eters and computational costs 0{N)). Therefore, we 
take one physically plausible Q in this study. 

By using the fe-dependent parameters tp^^\k) and the 
Gutzwiller factor Vq (eq. (21)), our variational wave 
function can also represent the resonating valence bond 
(RVB) basis, which is known to offer highly accurate 
variational wave functions in spin systems. We note that 
the RVB basis can represent the state with spin correla- 
tions decaying with arbitrary power laws for increasing 
distance. The relation between VQ\(ppa.iT) and the RVB 
basis is discussed in Appendix A. 

In quantum chemistry, Casula et al. have introduced a 
similar wave function called an antisymmetrized geminal 
power. Since the singlet pairs are only included in this 
wave function, it is not connected to the AF mean-field 
wave function. Our extension offers a clear representation 
to include the singlet pairing wave functions and the AF 
mean- field wave functions. 

For actual numerical calculations, we rewrite |(/>pair) in 
a real space representation: 



l^^'pair) 



E 



N/2 



|0) 



(18) 



with 



feeBZ 



(19) 



fceAFBZ 



Here, one of the parameters {(p^^\k), ip^'^\k)} is not in- 
dependent because of the normalization of the wave func- 
tion. 

2.1.2 Gutzwiller- J astrow factors 

In the variational study, the Gutzwiller- Jastrow type 
wave functions^^'^^ are often used to take account of 
many-body correlations. The Gutzwiller- Jastrow corre- 
lation factor V is operated to the one-body wave func- 
tion 1^), namely as "Pl^). Since \cp) is usually represented 
in the fe-spacc configuration, the factor V, constructed 
with many-body operator in the real space configura- 
tion, introduces compromise of real space and fe-space 
representations into one wave function. Because of this 
uncommutable nature, this factor allows us to go beyond 
the variational framework of a single Slater determinant 
and a linear combination of many Slater determinants 
are generated after the operation of V, which is crucial 
in representing strong correlation effects. In this paper, 
we adopt three many-body operators Vg, ^d4i' ^^'^ ^J' 



which are called the Gutzwiller factor, the doublon-holon 
correlation factor, and the Jastrow factor, respectively. 

Gutzwiller has introduced a basic and efficient corre- 
lation factor Vg,^^ which gives different weights to the 
wave function depending on the rate of double occu- 
pancy: 



Vg = exp 



-9 E "iT'^n = n ~ ~ ^ ^)"iT"a 



(20) 

where is a variational parameter. In the limit g 
oo, Vg fully projects out the configurations with finite 
double occupancy as 



=n[ 



1 - nifTiii 



(21) 



Vq is used for the Heisenberg model and the t-J model. 
In the Hubbard model with finite U/t, the double occu- 
pancy is nonzero even in the insulating state. Thus we 
deal with Vg at finite g. 

In order to take account of many-body effects beyond 
the Gutzwiller factor, the doublon-holon correlation fac- 
implemented in the wave function. This factor 
comes from the idea that a doublon (doubly occupied 
site) and a holon (empty site) are bound in the insulator 
for large U/t.^ The short-ranged correlation factor with 
many-body operators has a form 



Vd-h = exp 



"2E4'o) 



(22) 



where ai and are variational parameters. Here, 

is a many-body operator which is diagonal in the real 
space representations and may be given by 

"if a doublon (holon) exists at" 
the site i and m holons (dou- 
blons) surround at the ^-th . (23) 
.nearest neighbor 
f otherwise 1 



AO 

^z(m) 



1 



For example, is written by 

n.n. n.n. 
4o) = 11(1 - hi+r) + hi - di+r), 



(24) 



where the product H" " ^^^^ nearest-neighbor sites, 
and di = nifHii and hi = (1 — rii^){l — riii) are doublon 
and holon operators, respectively. The doublon-holon 
correlation factor Vd-h given by eq. (22) or by slightly 
different forms has been adopted in several VMC stud- 
Recently, there is a proposal to extend eq. (22) 



les 



36-40 



by introducing many ^l^^^y*^ We take V^\ as 



T'^i = exp 



where a|^j are variational parameters. It is in principle 
possible to include operators with m = 3, 4, . . ., but con- 
tributions of higher m parts have turned out to be negli- 
gible while have induced instabilities in our optimization 
procedure. Therefore we confine ourselves to m up to 2. 



EE"£)Ee 



m=0i=l,2 



(25) 
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Fig. 1. Weights of the many-body operators in Vd-h (upper row) 
and Vj (lower row) for some configurations. Filled circles, open 

circles, and filled diamonds denote doublon sites, holon sites, and 
single occupied sites, respectively. 



Jastrow has introduced a long-ranged correlation fac- 
tor for continuum systems. This factor takes into ac- 
count correlation effects through two-body operators. In 
the Hubbard model at quarter filling, Yokoyama and 
Shiba have discussed the effects of the Jastrow-type cor- 
relation factor. '^^ Recently, Capello et al. have claimed 
a necessity of this factor to describe the Mott transi- 
tion. ^°The Jastrow factor Vj in lattice models has the 
following form: 



Vj = exp 



(26) 



with two-body terms, where rii = rii^ + is a density 
operator and = u(rj — rj) are variational parameters 
depending on the displacement — Vj . The on-site Jas- 
trow factor is equivalent to the Gutzwiller factor except 
for a constant factor: 



= exp 



= exp 



-v{0)N 



exp 



const. 



(27) 

Prom the viewpoint of doublon-holon correlations, Vj 
can be rewritten as doublon-doublon (holon-holon) re- 
pulsive and doublon-holon attractive operators: 

VijUiUj = ^ Vij {didj + hihj - dihj - hidj) + C, 



■Pj oc exp 



^ Vij {didj + hihj - dihj - hidj) 



(28) 



where C is some constant. We remark the difference be- 
tween the doublon-holon correlation factor and the Jas- 
trow factor: For example, the two many-body operators 
give different weights to the configurations as shown in 
Fig. 1. 

2.1.3 Quantum-number projection 

In general, quantum many-body systems have several 
symmetries related to the Hamiltonian such as transla- 
tional symmetry, point group symmetry of lattice, U{1) 
gauge symmetry, and SU (2) spin-rotational symmetry. 
While symmetry breaking occurs in the thermodynamic 
limit, these symmetries must be preserved in finite many- 
body systems. 



Variational wave functions constructed from one-body 
parts and the Gutzwiller-Jastrow factors do not of- 
ten satisfy inherent symmetry properties, because the 
Hartree-Fock-Bogoliubov type one-body part comes from 
symmetry broken mean- field treatment. Even in the gen- 
eralized pairing wave function \(j)pair), the spin-rotational 
symmetry is broken by the orbital (/j^^^ (fe) which enables 
to include the mean-field AF state. 

The quantum-number projection technique^^ enables 
to control symmetries of wave function. This technique 
has been used successfully in the PIRG method^^ and the 
GBMC method. By using the quantum-number pro- 
jection together with the Gutzwiller-Jastrow factor, one 
can construct variational wave functions with controlled 
symmetries and many-body correlations. The quantum- 
number projection operator £ is constructed by super- 
posing transformation operators T^"^ with weights Wn'- 



(29) 



where \(f>) and |(/)^"^) are the original one-body part and 
the transformed one-body parts, respectively. When C 
restores some continuous symmetry, the summation J2n 
is replaced by the integration over some continuous vari- 
able. 

The SU{2) spin-rotational symmetry is restored by su- 
perposing wave functions rotated in the spin space. The 
spin projection operator jC^ which filters out = 
component of \^) and generates a state with total spin S 
and = has a form 

2S + 1 f 

= j dQPsicosP)R{Q), (30) 

where ^2 = (a, /?, 7) is the Euler angle and the integration 
is performed over whole range of Q. The weight P5 (cos (3) 
is the S'-th Legendre polynomial. The rotational operator 

R{Q) is defined as 



i?(/2) = R'ia)Ry{P)R''{"f) 



JaS' ^ifiS-" ^i-lS' 



(31) 



where 5^ and S'^ are total spin operators of y and z 
directions, respectively. 

Now we consider operating to the one-body part 
|(/)) which has the form 

N/2 



|0). (32) 



E E f^r 

.i,j=l ct,<t'=T,J. 

The rotated wave function R{f2)\(f>) is represented by the 
same form as eq. (32) with rotated creation operator 

elm 

1 



E E /r'4(^)sV(^2) 

U,j=l a,o-'=T,J. 



|0). 



(33) 

The rotated creation operator ct^(i7) is quantized along 
an ax;is rotated from z direction: 




R^(a)R^(/J)R^(7) 




(34) 
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je/2 









g-'ie/2 



cos(6'/2) sin(6l/2) \ 

-sin(6i/2) cos(6i/2) ) 

Thus, C^\(t>) has a form 
25+1 



(35) 
(36) 



dQPs{cos(3) 



N/2 



25+1 
87r2 



^ Kf'cimcl,m |0) 

_i,j,(T.(T 

da i dp I d-y sin P Ps {cos P) 















.lj,cr,<T 



JV/2 



|0), (37) 



where /^'^'^ (i?) is transformed from f[j°' by using eqs. 
(34), (35), and (36). 

The one-body part in this study introduced in 
§2.1.1 contains only 5^ = component |0) = 
[X^ij /yC-|C^J^/^|0). Then, the integration over 7 can 
be omitted and C^ldi) is written as 



25+1 
47r 



I da dp sin pPs {cos P) 
Jo 

X i?^(a)i?2'(/3)|</.). (38) 



In order to omit the integration over a, we rewrite cq. 
(38) to a more convenient form for VMC calculations. 
VMC is performed by sampling of a complete set of 
real space configurations {\x)} with which \ip) can be 
expanded: 



(39) 



Since the integration over a filters out 5^ ^ compo- 
nent, we can chose {\x)} with 5^ = condition and omit 
this integration: 



5^=0 

El 

X 

E 



25+1 



X 



dp sin pPs {cos p){x\Ry{P)\(l)). 

(40) 

When the Gutzwiller-Jastrow factor is operated to 
C^\(j)), a similar formula can be obtained as 



V£'\^)^ \x){x\V£' 



(41) 



with 

{xlVC'^l 



29+1 

P{x) — ^ I dp sin pPs {cos P) 

X {x\Ry{PM. (42) 
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Fig. 2. (Color online) Total spin {SP") as a function of the system 
size L at i'/t = 0, = 4, n = 1. Error bars are comparable to 
the symbol size. 
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Fig. 3. (Color online) Each element of the spin structure factor 
obtained by (a) PjP^rh^ol^pair) and (b) VjV^\VGCS^°\<t>p^ir) 
for t'/t = 0, U/t = 4, n = 1, and L = 14. Error bars are 
comparable to the symbol size. In (b), all the three cases are on 
the same curve. 



The integration over P is evaluated efhciently by the 
Gauss-Legendre quadrature in actual numerical calcu- 
lations. ''^ Typically, for 5 = of the half- filled electron 
system in L = 4 and L = 14 lattices, we need 10 and 20 
mesh points, respectively. 

Figure 2 shows the actual VMC simulation results on 
the total spin (5^) defined as 



(52) = ^(5..5, 



(43) 



The total spin (5^) grows as the system size increases if 
the projection is not imposed. On the other hand, spin 
projected wave function strictly keeps the 5 = state, of 
course. Figure 3 shows each element of the spin structure 
factor 

^"(9) = i-^(5r5;)e*'^-('^-'-^) {a = x,y,z), (44) 



S{q)^l{s-{q) + Sy{q) + S^{q) 



(45) 



The projected wave function recovers the symmetric 
property of 5"(q). 

Symmetry breaking of the spin part in |(/)pair) is caused 
by the components of y'^^-' (fc) , because ip'^'^\k) makes sin- 
glet and triplet pairs in |(/)pair)- There are two ways to 
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Table I. Variational energy of wave function with (a) only singlet 
pairs and (b) symmetry breaking orbitals and spin projection for 
L = 4, t' It = 0, U/t = 5, n = 1. The Gutzwiller-Jastrow factor 
is V = V^V^^Vq. The number in parentheses is the statistical 
error in the last digits. 



have 

{x\ct>) = {x\Y,{N/2)\ 



N/2 





Energy 


X 


{a)P|0pair) {¥'(^Hfe) = 0) 


-12.010(8) 




(b) Pr^=0|<^pai.) 


-12.459(6) 




Exa<;t diagonalization 


-12.5300 





^r-pf^2e-i)r-p(2t) ~ '"•p(2<)'"-p(2<- 



restore the spin rotational symmetry. The first way is 
restricting the pair orbitals only to ip'-^^k) and setting 
ip^'^\k) = 0. Then Impair) has only singlet pairs and the 
spin projection can be omitted. The second way is deal- 
ing with all orbitals (p^^\k), if^^\k) which requires an 
additional spin projection. 

Table I shows the variational energy calculated by the 
above two ways. The wave function on the latter way is 
better than one on the former. The combination of spin 
projection and symmetry breaking pair orbitals provides 
much accurate wave functions. 

There are other quantum-number projections to re- 
store symmetries. The total momentum projection and 
the lattice symmetry projection restore the translational 
symmetry and the point group symmetry of lattice, re- 
spectively. These projections arc not adopted in this 
study. The reason is the following: Although the one- 
body part in cq. (16) can break the translational sym- 
metry (formation of A-B sublattice) by the AF symme- 
try breaking, this broken symmetry can be restored by 
performing the spin projection, because a superposition 
of the rotated wave functions includes a superposition of 
wave functions translated between A and B sublattice. 
Therefore, after the spin projection, we do not need to 
introduce other types of symmetry restoration. 

2.2 Calculation of inner product {x\ip) 

For actual VMC calculations, the inner product be- 
tween a real space configuration |.x) and a given wave 
function \ip) is a key quantity. In this section, we explain 
the way to calculate {x\ip). First, we explain the inner 
product between \x) and one-body part \(f>) in a general 
case with N particle {N: even). The real space configu- 
ration I a;) has the form 

4 



4 



.10) 



(46) 



and one-body part \<j)) has a general pairing functional 

form 



E 

.i,i=l <T,cr' = T.i 



ct ct , 



N/2 



|0). (47) 



Here, Ff^"' is, for example, given by fij defined in eq. 
(19). As first pointed out by Bouchaud et al.,^^ the in- 
ner product {x\(l)) is given as a Pfaffian oi N x N skew- 
symmetric matrix. We explain this fact below in detail. 
Expanding |0) and picking up nonvanishing terms, we 



n( 

e=i 

N/2 

n (""»'T'(2«-i 



)(7-p(2e-l) ^V(2i)<^V{2t) 



-.;■) 

)|0), (48) 



where V is the permutation of A'' indices with the condi- 
tion 



r{2i-l) < V{2i) 

■P(l) <-p(3) <■■■ <V{N -I) 



(49) 



From the power N/2 in the term with the same 
element of Y{{F^;,' - F^:^^)Y{{cl^cl,^,) appears {N/2)\ 
times. The commutation relation of fermion operators 
gives the sign 

N/2 

n ('^'■P(2<:-1)<^T'(2*-1)^'-P(2*)<^P(2*)) = ' (50) 

where (—1)^ is the parity of V. Thus, (a;|(/)) has the form 
(x|</.) = (Ar/2)!^(-l)^ 



N/2 

£=1 



pCr-p(2«-l)0'P(2f) 
'■P(2f-l)'"p(2f) 



pfp(2f)fp(2f 



-.;■) 



(A^/2)!PfX, 



(51) 



where Pf X is a Pfaffian oi NxN skew-symmetric matrix 
X with the element 



The linear algebra of skew-symmetric matrix and Pfaf- 
fian is described in Appendix B. 

Next, we consider the inner product of a spin pro- 
jected wave function with the Gutzwiller-Jastrow factor 
V: = VC^\<j)). From eq. (42), the inner product is 
given by 

29-1-1 

{x\PC^\(t)) = P{x) — / dp sin /3Ps (cos /3) 
2 Jo 

X {x\RyiPM, (53) 
where the condition = is imposed to \x): 
la;) = i,c^ * • • • rC^ i---c^ 1 10), (54) 

1-^/ ''riT r2T '-JV/2T I'N/2+li f„/2+2i fivil"/' V" 

and |(^) has a form 

-1 N/2 



|0). 



(55) 



8 J. Phys. Soc. Jpn. 



Full Paper 



D. Tahara and M. Imada 



The component (a;|i?^(/?)|(?i) can be evaluated by 

{x\Rym<t>) = {A 



are the components of the energy gradient vector g, and 



^/,,(cos(/?/2)cl^+sin(/3/2)clJ 
X (-sin(/3/2)ct^+cos(/3/2)c]^) 



N/2 

|0) 



{x\ 



-1 N/2 



|0) 



= PfX(/?), (56) 

where the element of skew-symmetric matrix X(/3) is 

f F}},^ iP) - FH , (/?) (. < iV/2, J < N/2) 

(/3) - F}J,^ (/?) (z < 7V/2, J > iV/2) 

F4^. (/?) - {13) (z > N/2, J < N/2) " 

[ , - FJ.i^ (z > 7V/2, j > 7V/2) 

(57) 

Therefore, the inner product {x\VjC^\(p) is obtained as 



P{x) 



2S+1 



dp sin /3Ps (cos /3)PfX(/3). 

(58) 



3. Optimization Method 

In the procedure of VMC, the wave function optimiza- 
tion is one of the most important tasks. In the optimiza- 
tion, we have to keep in mind the following limitations. 

(i) The estimated value of the cost function (usually 
the total energy) and its derivatives have the sta- 
tistical noises by MC samplings. 

(ii) There is a trade off between computational costs 
and accuracy when one employs the estimation of 
higher-order derivatives of the energy in the varia- 
tional parameter space. 

In this chapter, first we summarize the basic idea of wave 
function optimizations by energy minimization. Then, 
the stochastic reconfiguration (SR) method, ^'^ which 
Sorella has developed in order to optimize many param- 
eters, is explained in detail. 

3. 1 Basic idea of wave function optimization — Steepest 
Descent method and Newton method 

We discuss an efficient way of minimizing the en- 
ergy Ea = (ipal'Hltpce) / {tpali^a} estimated from the 
wave function \ipct) with variational parameters {ak\k = 
I,-- - ,p}. Here a denotes the initial vector in the p- 
dimensional parameter space. 

The energy Ea+-y is expanded up to the second order 
around a: 

p 1 ^ 

E„,+^=Ea, + Y,9k7k + ^ hulkli + Oi-i^), (59) 

fe=l k,l=l 

where 7 is the vector for parameter variations, 
d 



9k 



dak 



-En 



(fc = l,--- ,p) 



(60) 



hu = 



92 



-E^ {k,£=l,--- ,p) (61) 



dakdai 

are the elements of the energy Hessian matrix h. 

With the first order approximation, the steepest de- 
cent (SD) method gives the updated variational param- 
eter by 

a'k = ak+Jk, (62) 
where the change from the initial value ak should be 



Ik 



-Atgk {^ = -Atg). 



(63) 



Here, At is a small constant. Combination with the 
second-order information, i.e. the Hessian, leads to the 

Newton method. By imposing the stationary condition 
dEa/dak = {k = 1, - ■ ■ ,p), the best parameter change 
7 is obtained by 

p 

lk = -Y.Kht (7 = -h-^9)- (64) 
£=1 

Let us generalize eqs. (63) and (64) with suitably cho- 
sen nonsingular matrix X: 



Ik 



p 



Xki'gt (7 = -x-^g). 



(65) 



Equations (63) and (64) are reduced from eq. (65) by 
setting Xki = Ski/ At and Xke = hki, respectively. 

As long as the energy gradient g is estimated with the 
mathematically correct formula, the parameter change 
7 will converge at the minimum or at a stationary point 
irrespective of the choice of X. However, the computa- 
tional efficiency strongly depends on the choice of X and 
X should be chosen to accelerate the optimization within 
computational stability. 

3.2 Stochastic Reconfiguration method 

Sorella has developed the SR method, which offers a 
simple but very stable optimization method. In order to 
deal with a large number of variational parameters and 
optimize all the parameters simultaneously, we employ 
the SR method in this studies. First, we introduce the 
normalized wave function 



1 



(66) 



Then the expansion of |V'a+-y) up to the first order 
around a is 

_ _ ^ _ 

|V^«+7) = + + 0(7'), (67) 

k=l 

where IV'fea) {k = - ■ ■ iP) are the derivatives of |V'a): 
d - 

l-t/jka) = ^ IV'a) 

oak 

(68) 
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The wave function set {|V'fca)|^ = I,''' ,p} forms 
nonorthogonal basis in the p-dimensional parameter 
space. The norm of the variation between |^c«) and 
IV'a+-r) is 

2 



= X] IkJeii'kali'ea) = X! ^kleSke- (69) 



k,e=i 



kj=l 



Since Ski = {'>Pka\4'eoc) is the overlap matrix in the pa- 
rameter space, S becomes positive definite even with a 
finite number of samples. Equation (69) shows that S is 
the metric matrix in the parameter space. 

The SR method chooses S as the matrix X in eq. (65), 
namely 



7fc 



-Atj2^ki9i {l = -AtS-'g), (70) 



where At is a small constant. As the overlap matrix S 
does not have any information about the energy Hessian, 
the SR method is close to the SD method. The main dif- 
ference is that the SR method takes into account the 
variation of the wave function. We can derive eq. (70) 
by minimizing the functional !Fs^ = AE^^^ + XA'^^^.^ 
with a Lagrange multiplier A. Here AEn^. = X^j. Qk^k is 
the linear change of the energy. The stationary condition 
dJ^SK/dlk = {k = 1, • • • ,p) leads to the SR formula 
(70) with At = (2A)-^ The SD method can be obtained 
in a similar way. We can derive eq. (63) with At = (2A)~^ 
by minimizing the functional J-su = AEnn. + AzigQ with 
ZigQ = X]fe7fc, where Asu is the Cartesian distance in 
the parameter space. The advantage of the SR method 
compared with the SD method is the following: Some- 
times small change of the variational parameters cor- 
responds to a large change of the wave function, and 
conversely a large change of the variational parameters 
corresponds to a small change of the wave function. This 
leads to uncontrolled changes of the wave function if one 
takes Ske = Skt (SD method). When the change of the 
wave function exceeds a threshold, the iteration for the 
optimization becomes unstable. To suppress this insta- 
bility, one needs to keep At small enough for the event 
of the largest change of the wave function. If one can con- 
trol change of the wave function, At can be taken large. 
The SR method takes into account this effect through 
a better definition of the distance /\norm- Thus, the SR 
method is more stable than the SD method. Finally, we 
can choose larger At to accelerate the convergence. 

3.3 Stabilization of SR method 

In the VMC calculation, it is important to optimize 
the wave function stably with a small number of samples. 
The main instability in the SR method comes from the 
overlap matrix S. Though the overlap matrix S is positive 
definite even with a finite number of samples, the inverse 
matrix amplifies the statistical noise in the energy 
gradient g when the ratio of the maximum eigenvalue 
and the minimum eigenvalue becomes extremely large. 
The statistical noise of MC sampling and a variety in 
dependence on each parameters enlarge this ratio. 



To stabilize the SR method, we apply two tech- 
niques:^^' 

(i) Modification of diagonal elements in S, 

(ii) Truncation of redundant directions in the parame- 
ter space. 

3.3.1 Modification of diagonal elements in S 

As we have already mentioned, deformation of X in eq. 
(65) does not change the optimal parameters. Then we 
follow the stabilization method in ref. 44 by modifying 
diagonal elements in S: 

Skk^{l + e)Skk, (71) 

where £ <C 1 is a small constant. This modification pre- 
serves the positive definite property, because the sum of 
two positive definite matrices Ske and sdkeSke remains a 
positive definite matrix. This modified matrix pulls up 
extremely small eigenvalues and suppresses fluctuations 
in the SR iteration. 

As in the stabilization of the ordinary Newton method, 
it is possible to add a uniform constant e to diagonal 
elements {Skk Skk + £)• This also stabilizes the SR 
method, but the convergence becomes slower than the 
former, because this modification does not take account 
of the metric in the parameter space. Though this stabi- 
lization in eq. (71) suppresses fluctuations in a large part 
of the variational parameter space, it becomes inefficient 
in some cases. When the ratio of the maximum value and 
the minimum value of the diagonal elements becomes ex- 
tremely large, the ill part is not modified efficiently by 
the additional matrix sSkeSke. In order to stabilize the SR 
method with large statistical noises further, we combine 
the modification in eq. (71) with a truncation technique 
as discussed in the next part. 

3.3.2 Truncation of redundant directions 

Casula et al. have introduced a truncation technique 
for irrelevant variational parameters to stabilize the SR 
method. '^^ They have directly truncated some parame- 
ters. Here, we introduce a better truncation procedure 
by considering the eigenvalues of the overlap matrix S. 

As S is a p-dimensional positive definite symmetric 
matrix, we can diagonalize S by an orthogonal matrix U: 



p 

E 

k,e=i 



Ski — ^ KUkiUii 



•^ke. 



p 



1 



X -T-UkiUii 

1=1 



(72) 

where A^ > (i = 1, • • • ,]?) are the eigenvalues and ar- 
ranged in descending order (Ai > A2 > • • • > \p). Then 
we apply an orthogonal transformation to the parameter 
change 7: 

p p 
Xk = ^1kUki Ik = ^UkiXi. (73) 

k=l i=l 

From eqs. (69), (72), and (73), we obtain the following: 



^norm = E IkliXtUkrUi, = 



(74) 



k,e,; 



i=l 
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This means that the variation in the direction Xi where 
Ai/Ai < £wf is satisfied is redundant in the range of rela- 
tive accuracy £wf- Moreover, this direction brings an in- 
stability (~ l/£wf) to S~^. Therefore, we can control the 
relative accuracy of ^noim ^^'i stability of S^^ by tuning 
Evff. We truncate the direction {xi\i = q+1, - ■ ■ ,p} which 
satisfy Aj/Ai < Ewf- The SR formula (70) is changed into 

p p 



7fc = -At S^lgt = -^tYY^ Y^kiUuge (75) 

fc=i i=i 
J| (Truncation) 



1 

Ik = -Aty^^—UkiUeigi. 



(76) 



=1 i=l 



By introducing the truncation parameter e„f, wc can 
control both the accuracy and the stability of the SR 
optimization. The truncation by decomposing S into the 
orthogonal matrix U (eq. (72)) is essentially equivalent 
to the singular value decomposition (SVD) method''^ to 
handle ill matrices. In the optimization method of vari- 
ational wave functions, the SVD has been adopted in 
the linear method developed by Nightingale and Melik- 
Alaverdian, which uses some information of the Hes- 
sian.25.27 

3.3.3 Practical parameters in SR method 

In order to perform the stable SR optimization with 
small MC samples (4 x 10^-10'* samples), we practically 
choose the parameters in the SR method as At = 0.1, 
e = 0.2, and Swf = 0.001. The optimized variational pa- 
rameters are typically obtained by averaging the param- 
eters over 100 optimization steps after 500-5000 steps. 

3.4 Derivative operator 

In order to perform the SR optimization, we need 
to evaluate the overlap matrix S and the energy gra- 
dient g. We follow a standard way in continuum systems 
by introducing derivative operators. ^^'^^ In lattice mod- 
els, Sorella have derived representations of these opera- 

tors.23.24 

The wave function considered here is a quantum- 
number projected one with the Gutzwiller-Jastrow factor 

\^c)=Vc.C\M=VcJ2'^n\^^^^), (77) 

n 

where the summation 5^,^ and weights w„ come from 
the quantum-number projection C. Here, the Gutzwiller- 
Jastrow factor Va and the transformed one-body parts 
\(j)a^) have variational parameters. We assume that the 
parameter dependence in Va and \(l)a^) is separated. 

First, we introduce operators Ok (fc = 1, • • • ,p) which 
are diagonal in real space configurations \x). is defined 
as 



El-) 



1 



d 



{x\ipa) dak 



{x\lpcx) 



x)Ok{x){x\ 
(78) 



and satisfies the following relations: 

{x\Ok\i^c.) = ^(a^lV-c), 

The derivative of the normalized wave function (eq. (68)) 
is rewritten as 



(79) 
(80) 

(81) 



-.{Ok-{Ok)Moc), 



(82) 
(83) 



and the overlap matrix S is obtained as 

Ske = {OkOe) - {Ok){Oe). 
The energy gradient g is also written as 

gk = -^(V'c«|W|V'«) = (V'fcc.lWIV'a) + (V'clWIV'fea) 

dak 

= 2{HOk)-2{H){Ok). (84) 

Next, we derive expressions for Ok{x). The Gutzwiller- 
Jastrow factor has an exponential form 



Vex = exp 



(85) 



where ak are variational parameters and 0k are diag- 
onal operators for real space configurations {0k\x) = 
6'fe(a;)|a;)). Ok{x) corresponding to Va is obtained as 

Ok{x) = -Ok{x). (86) 

Since the inner product {x\(pa^) is evaluated by the Pfaf- 
fian Pf X^"*, Ok{x) corresponding to the one-body part 
is related to the derivative of the Pfaffian. The result is 
given as 

1 a , , , , 1 d 



Okix) 



{x\ipa) dak 
1 d 
{x\C\(j)a) dak 

1 



{x\C\(l)a} dak 



{x\C\(j)a) 



E 



{x\C\(t)a) 



^""dak " 



" dak " 



n) 



(87) 



The overlap matrix S and the energy gradient g are eval- 
uated under the VMC sampling by using the above for- 
mulae. For example, if \(j)a) is given by |(^pair) defined by 
eq. (16), dXa^/dak contains derivatives dfij/dip^^\k) 
and dfij/dip'^'^\k) through eq. (19). 

4. Results 

In this section, we apply our improved wave functions 
to variational calculations for the two-dimensional Hub- 
bard model. We compare variational results obtained 
from our wave functions with unbiased results obtained 
by the ED method and the AFQMC method.'^ 
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Table II. Comparison of variational energies obtained from different wave functions and the exact energy calculated by ED for L = 4, 
t'/t = 0, and n = 1. The numbers in parentheses are the statistical errors in the last digits. 





U/t = 4 


5 


10 


20 


35 


I0af; 


— iz.UZo[\J ) 


— lU.y ( y t^Uj 


— D.lUoy(^UJ 


— o. iOiUl^Uj 


— i.oZlZiiJ ) 


' G 1 <PAF / 


1 A no/i ^ 
— L^.VJyl I 


1 1 ao/'o^ 

— Ll.oZyZ ) 


a /1 1 /I ^ 
— 0.41(^1 ) 


— o.z^yL ) 


— l.o^yz ) 


7'd-h7'Gl0AF> 


— 14.204(1) 


— 11.951(1) 


—6.551(1) 


—3.42(1) 


— 1.978(5) 


.CS=O|0AF> 


-13.735(5) 


-11.775(5) 


-6.72(1) 


-3.54(1) 


-2.04(2) 


7'd-hPG£^="l0AF> 


-14.437(2) 


-12.356(1) 


-6.946(1) 


-3.63(1) 


-2.10(1) 


^J^d"h^G.CS=O|0AF> 


-14.443(2) 


-12.363(2) 


-6.952(4) 


-3.629(6) 


-2.101(4) 


^J^d^h^Gl-jipair) 


-14.278(1) 


-12.10(1) 


-6.738(7) 


-3.486(7) 


-2.009(8) 




-14.512(3) 


-12.459(6) 


-7.050(2) 


-3.693(2) 


-2.135(4) 


Exact (ED) 


-14.5935 


-12.5300 


-7.13239 


-3.76124 


-2.18092 



Table III. Energy E, double occupancy D, and spin structure factor S{n, n) of ED and VMC results. The variational wave function is 
^J^d?h^G'C'^~''l'/'pair)- The numbers in parentheses are the statistical errors in the last digits. 





t'/t = 






t'/t = -0.3 






t'/t = -0.5 






U/t 


4 


8 


10 


4 


8 


10 


4 


8 


10 


E-ED 


-14.5935 


-8.63871 


-7.13239 


-14.6322 


-8.66093 


-7.15550 


-14.7305 


-8.71446 


-7.21193 


EvMC 


-14.512(3) 


-8.569(2) 


-7.050(2) 


-14.531(2) 


-8.571(2) 


-7.048(3) 


-14.561(3) 


-8.578(2) 


-7.017(6) 


Deo 


0.14389 


0.05664 


0.03911 


0.14410 


0.05661 


0.03910 


0.14463 


0.05650 


0.03906 


DyMC 


0.1441(3) 


0.0570(1) 


0.0386(2) 


0.1442(2) 


0.0569(3) 


0.0394(4) 


0.1448(3) 


0.0573(2) 


0.0413(2) 


SedC'^i't) 


0.63517 


1.21722 


1.31045 


0.61656 


1.20136 


1.28914 


0.56647 


1.15483 


1.22434 




) 0.674(3) 


1.250(2) 


1.368(2) 


0.664(1) 


1.260(2) 


1.367(3) 


0.648(1) 


1.261(2) 


1.371(4) 



0.20 1— n — r 




Fig. 4. (Color online) Relative accuracy of different wave func- 
tions for L = 4, t'/t = 0, and n = 1. The accuracy zig = 
1 — SvMc/^ED denotes relative difference of the variational en- 
ergy -EvMC the exact value Bed- 



Comparison with exact diagnalization and energy 
variance 

First, we compare the energy of different variational 
wave functions with the ED results. Table II shows the 
results for the system L = 4, t'/t = 0, U/t = 4, and 
n = 1. The one-body part denoted by |^af) is the mean- 
field state which diagonalizes the Hamiltonian (5) with 
^scC^) ~ ^sc(^) ~ ^^'^ with ZiAF being optimized 
as a variational parameter of |^af)- 

Figure 4 shows the relative accuracy of the above re- 
sults. The spin-projected generalized pairing wave func- 
tion with the Gutzwiller-Jastrow factor has the best ac- 
curacy. The spin projection acts efficiently in each wave 
function. Restoration of spin rotational symmetry and 
filtering out of excited states with other spin quantum- 



Table IV. Comparison of energies and energy variances obtained 
from wave functions with different components and the exact 
energy calculated by ED for L = i, t' /t = 0, U/t = 4, and 
n = 1. The numbers in parentheses are the statistical errors in 
the last digits. 





E 


{{■H') - {W)/{nr 


Impair) 


-13.006(6) 


0.0350(2) 


.C^ = °l</'pair> 


-13.831(4) 


0.0231(6) 




-14.278(1) 


0.0096(3) 




-14.489(2) 


0.0032(2) 


•Pd?h^G/:^ = °l</'pair) 


-14.509(2) 


0.0028(1) 


^j7'|rh^G^^ = °l</'pair> 


-14.512(3) 


0.0028(3) 


Exact (ED) 


-14.5935 






numbers are crucial to improve variational wave func- 
tions. The Jastrow factor and the extension of the 
doublon-holon correlation factor do not seem to lower 
the energy substantially. Table III shows the compari- 
son at nonzero t'/t for the same lattice size. The accu- 
racy slightly declines as the frustration increases. The 
double occupancy D in Table III is defined as D = 

(i/7Vs)E.K"a)- 

Table IV shows the comparison of the energy and the 
energy variance {{'U?) - (H)'^) / {H}'^ . It shows that the 
improvement is brought step by step by implementing 
each refined treatment of Impair), >C'^"°, Vg, 'Pd^hi and 

4-2 Size dependence 

Next, we examine size dependence by comparing with 
AFQMC results, which are exact within the statistical 
accuracy. Figure 5 shows the size dependence of total 
energy per site E/Ng, whose difference from the ther- 
modynamic limit L oo scales proportional to as 
is derived from spin-wave theory."*^ The spin-projected 
generalized pairing wave function with the Gutzwiller- 
Jastrow factor has the best variational energy. The en- 
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Fig. 5. (Color online) Total energy E/Ns as a function of 

for t' /t = 0, U /t = 4, ?i = 1. The exact values are calculated by 
ED (L = 4) and AFQMC (L = 6,8,10,12,14). Error bars are 
comparable to the symbol size. 



ergy gain by the spin projection becomes smaller when 
the system size increases, because the spin projection 
lowers the energy by filtering out higher-energy excited 
states with other spin quantum-numbers while the en- 
ergy is not efficiently lowered when the energy of exci- 
tations belonging to the same quantum number as the 
ground state becomes close to the ground-state energy 
in larger sizes. In finite size calculations, however, the 
spin projection is still useful to obtain better wave func- 
tions on a quantitative level. Furthermore, the energy 
of the spin-projected wave functions is better scaled by 
than that of unprojected cases, which makes the 
extrapolation to the thermodynamic limit easier. 

Introducing a large number of variational parameters 
into the one-body part efficiently improves the wave func- 
tion and allows us to go beyond the simple AF order 
generated from the mean-field treatment. We remark 
that the additional parameters for |0af) 



Ag^^\coskx — cos ky) and fi in eq. (5)) do not improve 
the wave function in this parameter region. We confirm 
this fact in the small system L = 4 and L = 6. 

Figure 6 shows the size dependence of the double oc- 
cupancy. |(/)pair) have smaller value than \<Paf)- This dif- 
ference mainly improves the variational energy. Figure 
7 shows the size dependence of spin structure factor 
S{tt,tt)/Ns plotted as a function of 1/L following the 
scaling by the spin-wave theory.^^ The behavior of spin- 
projected wave function is qualitatively better than that 
of the unprojected case. The staggered magnetization in 
the thermodynamic limit mg — [lim^Vj^^oo *S'(7r, 'k)/Ns]^^^ 
of each result is estimated as follows: 



(88) 



0.173 ± 0.004 : VMC |0af) 
0.158 ± 0.006 : VMC j^pah) 
0.138 ±0.005 : AFQMC^^ 
[ 0.304 ± 0.004 : Heisenberg model-*^) 

Though the variational results have a general tendency 
of showing larger S{Tr,Tr)/Ns and rus than the AFQMC 
results, the value obtained from |(/ipair) shows substantial 
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Fig. 6. (Color online) Size dependence of double occupancy. The 
exact values are calculated by the same way as Fig. 5. 
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Fig. 7. (Color online) Size dependence of spin structure factor. 
The exact values are calculated by the same way as Fig. 5. 



improvement as compared to that obtained from |^af)- 
This means that Impair) treats the AF correlation and 
quantum fluctuations more correctly than the mean-field 
descriptions, which often overestimate orders. 

From the above results, the variational wave function 
constructed from the combination of Impair), spin pro- 
jection, and the Gutzwiller-Jastrow factor offers the best 
description of the ground-state wave function among var- 
ious choices of variational functions. 



^.3 Spin correlation in hole-doped systems 

In order to examine accuracy under severe conditions, 
we calculate the spin correlation in the hole-doped sys- 
tems with t'/t ~ and U/t = 4. Figure 8 shows the 
doping dependence of the peak value of spin structure 
factor S'(qpcak). Although we are not able to compare our 
VMC results directly with the AFQMC results in ref. 7 
because of the difference in boundary conditions, our re- 
sults show excellent agreement with the unbiased results 
in a wide range of doping concentration 6 ~ 1 — 7i. This 
suggests that quantum fluctuations, especially the short- 
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Fig. 8. (Color online) Peak value of spin structure factor S'(cfpeak) 
as a function of the doping concentration S for t' /t = and U/t = 
4. Diamonds are AFQMC results reported in ref. 7. Triangles and 
circles are VMC results obtained by using the variational wave 
function 'Pj'P™h'^G'C^=°|0pair>- Error bars are comparable to 
the symbol size. The solid curve is the fitting given in Fig. 13(a) 
of ref. 7. The fitting satisfies S'(qfpeak) <^ for ^ ^ 0.2. 



ranged AF correlation which is cruciaUy important in 
hole-doped systems, are satisfactorily taken into account 
in our variational wave function. In fact, the system size 
dependence of S'(<7pcak) is small for S > 0.05 in agreement 
with AFQMC data.^ Our VMC and AFQMC data both 
consistently show that the AF long range order is re- 
stricted to the doping region much smaller than S = 0.05 
with the scaling 5(qpoak) oc S^^ in the paramagnetic re- 
gion (S > 0.05). 

4-.4 Other results 

We also apply our improved variational wave func- 
tion to the frustrated Hubbard model with t' /t = —0.3 
and n = l.'*^ The Mott transition between the para- 
magnetic metal and the AF insulator takes place at 
Uc/t — 3.3 ± 0.1 which can be favorably compared with 
the estimation by the PIRG method (Uc/t ~ 3.6). This 
Uc/t is much smaller than the previous variational esti- 
mate {Uc/t ^ 6.7).'*" The double occupancy keeps a large 
value {D 0.2) in the metallic phase near the Mott 
transition and shows very small U/t dependence. In the 
previous studies, this characteristic feature has been ob- 
tained only in the PIRG results. The variational wave 
functions employed in the literature include many-body 
correlations only by much restricted forms, such as the 
doublon-holon short-ranged factor. Such restricted form 
does not sufficiently take into account quantum fluctua- 
tions, which are strongly enhanced around the Mott tran- 
sition. Introducing a large number of variational param- 
eters in the Gutzwiller-Jastrow factor as well as in the 
one-body part allows quantitatively accurate treatment 
of fluctuations with complicated correlations. 

5. Summary and Discussions 

In this paper, we have extended the variational 
Monte Carlo (VMC) method and applied it to the two- 
dimensional Hubbard model. The VMC method origi- 
nally has several advantages for studies of the strongly 
correlated systems. This method is tractable in large 
system sizes even with strong interactions and geomet- 
rical frustrations. However, the bias inherently and in- 



evitably contained in the assumed variational form of the 
wave functions is a fundamental drawback in the VMC 
method. 

In order to overcome and go beyond the conventional 
limitation in the VMC framework, we have improved 
variational wave functions by the following extensions: 

(i) By introducing a large number of variational pa- 
rameters, we have constructed the one-body part 
including various states such as paramagnetic met- 
als, antiferromagnetically ordered states, and su- 
perconducting states with any wavenumber (spa- 
tial) dependence of gap functions within a single 
functional form. Moreover, Impair) enables efficient 
treatment of quantum fluctuations of spins. This 
extension efficiently reduces biases coming from as- 
sumed states in the previous VMC method. 

(ii) We have introduced a new factor for the variational 
wave functions: The quantum-number projection 
factor restores the inherent symmetry of the wave 
function and, as a result, the accuracy is substan- 
tially improved. 

(iii) We have combined our improvements with the re- 
cently improved Gutzwiller-Jastrow factor includ- 
ing many variational parameters. In particular, the 
fc-dependent Jastrow factor efficiently takes into ac- 
count fluctuations of charge. 

The accuracy of our variational framework has been 
examined by the comparison with the unbiased re- 
sults obtained from the exact diagonalization and the 
auxiliary-field quantum Monte Carlo method. It has 
turned out that the improvement of the one-body part 
and the quantum-number projection enable highly ac- 
curate descriptions of the wave function for the ground 
state. In the system with t' /t — 0, U/t — 4, and n — 1, 
the relative error reaches as low as 0.5% at L = 4 and 1% 
in the thermodynamic limit. These errors are typically a 
half of the best available results in the literature. 

Our improvement of the variational wave function does 
not change the basic numerical framework of the VMC 
method, namely sampling of the real space configura- 
tions. Therefore, our approach is able to combine the 
"post- VMC" method such as the Lanczos method^" and 
the diffusion Monte Carlo method.^* These additional 
treatments certainly even more reduce the biases. In this 
paper, we have used a single pairing wave function for 
the core one-body part. The linear combination of several 
pairing wave functions for the core will improve the vari- 
ational wave function, though this causes linear increase 
of the computational costs. The search for efficient repre- 
sentations of this multi-configuration is a future problem. 

Our improvement of variational Monte Carlo method 
opens a possibility of studying strongly correlated elec- 
tron systems by reducing the effects of biases from re- 
stricted variational forms. In particular, effects of short- 
ranged spin and charge fluctuations may be studied with 
quantitative accuracy. Applications of our refined algo- 
rithms will be reported elsewhere. 
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Appendix A: Relation between 7'^|</>pair) and 
the RVB basis 

The wave functions based on the resonating valence 
bond (RVB) basis'^'* provide highly accurate descriptions 
in spin systems. As discussed by Liang, Doucot, and An- 
derson,^'' these functions can control AF correlations. It 
is worth clarifying the relation between the RVB wave 
function and wave functions used in itinerant electron 
systems. As pointed out first by Anderson, the RVB 
wave function is equivalent to the so called "projected 
BCS" wave function. 

The RVB wave function is constructed by singlet- 
dimer covering as 

-7V/2 

n '^■"'3>' 

.k=l 



with 



10: 



RVB) - ; , 

{C} 



|0), (Al) 



where the summation X]{c} ''uns over all the patterns 
of coverings {C} with some conditions described below. 
When the AF order has A and B sublattices, two condi- 
tions "ifc e A" and "jfc e B" are imposed on {C}. Since 
two singlet dimers sharing one same site make a doublon: 

(* 7^ J> i ^ k ^ i), \(j)-RyB) can be rewritten as 

-| N/2 



t J 



N/2 



|0), 



|0) 



(A-2) 



where Vq 



J^j(l — ni-\nii) is the Gutzwiller factor (eq. 
(21)) which filters out all the components with a finite 
double occupancy, j) {c} sums up all the pairs 
appearing in {C}, and off-diagonal elements of are 
related to hi^ as 







(i € A, j € B) or {i € B, j € A) 

{i e A, j e A) or (?; e B, j e B)' 

(A-3) 

The diagonal elements fa are arbitrary because of Vq. 
Equation (A-3) means that /y = /(r^ — rj) depends on 
the relative vector (r^ — rj) and the translational sym- 
metry is preserved. Figure (A-1) shows three examples of 
A-B sublattice patterns and some relative vectors corre- 
sponding to nonvanishing elements /(t). By using the 
Fourier transformation, |0rvb) is written as 

-1 N/2 



10: 



RVB, 



E'^'^^T^- 



|0) = |p-BCS) (A.4) 



A = E/w 



-ik-r 



(A-5) 



Since [X^fe /fcCfc^cl^J^/^lO) is the BCS wave function, 
this means that the RVB wave function corresponds to 
the BCS wave function with Vq, which is called a "pro- 
jected BCS" wave function |p-BCS). 

From the above discussion, the variational wave func- 
tion constructed from |(/)pair) and Vq includes the RVB 
basis. Determination of fij or hij by hand is discussed 
in detail in the literature. ^'^ In our calculation, fij is nu- 
merically determined by the optimization of variational 
parameter ip^^^{k) in eq. (16). 












(b 


'(tt.O) 




a) {it, tt) 



Fig. A-1. (Color online) Possible A-B sublattice pattern of (a) 
(tt, 7r)-order, (b) (tt, 0)-order, and (c) (-n-, •n-/2)-order. Open circles 
and filled circles are A and B sublattice respectively. Vectors are 
some examples corresponding to nonvanishing elements /(r). 



Appendix B: Linear Algebra of Skew- 
symmetric Matrix and Pfaffian 

In this appendix, some useful relations between the 
skew-symmetric matrices and the Pfaffians are derived. 
We assume that all the elements of the matrices are real 
numbers. 

B. 1 Definition 

A 2M X 2M skew-symmetric matrix A satisfies the 
following relation: 

A^ = -A {A,j=-A,i), (B-1) 

where A-^ denotes the transposed matrix of A. The Pfaf- 
fian of A is defined as antisymmetrized product 

Pf A = ^[^12^34 • • • A2M-iaM] 
M 

= E(-i)^n^n2fe-im2fe)> (B-2) 

V k=l 

where the sum runs over all the pair partitions V of 2M 
indices such that 'P{2k - 1) < 7^(2^). Here, (-1)''' is the 
parity of the permutation V: 



1 2 
V{1) V{2) 



2M-1 2M \ 

V{2M~1) V{2M))- ^ ' 



B.2 Identities 

B.2.1 Basic relations 

The Pfafhan satisfies the following relations: 

Pf A'^ = (-l)*^Pf A (B-4) 
(PfA)2 = detA (B-5) 
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Pf 



Ai 
A2 



Pf Ai X Pf A2 



Pf (BAB'^) = detB X Pf A 



Pf 



C 

-C^O 



(-l)^(^-i)/2detC 



(B-6) 
(B-7) 
(B-8) 



where A, Ai, A2 are 2M x 2M skew-symmetric matrices, 
B(2M X 2M) and C(M x M) are arbitrary real matrices. 

By taking B as elementary transformation matrices, 
we can verify the following properties of PfafBans: 

(i) Multiplication of a row and corresponding column 
by a constant is equivalent to multiplication of orig- 
inal Pfaffian by the same constant. 

(ii) Interchange of two different rows and corresponding 
columns changes the sign of Pfafiian. 

(iii) Consider a vector p which is the same as a row of 
A and another transposed vector which is the 
same as the corresponding column of A. Addtion 
of p to another row and addition of to the cor- 
responding column does not change the value of 
PfafBan. 

By performing a Gaussian elimination technique and piv- 
oting rows and corresponding columns, wc can transform 
any skew-symmetric matrix into block-diagonal form and 
obtain the value of PfafBan: 



Pf A= (-l)Ppf 



Ai 
-Ai 



A2 
-A2 



Am 
-Am 



= (-If AiA2---A 



M, 



(B.9) 



where p is the frequency of pivoting and Ai {i = 
1, • • • , M) are results of the elimination. The calculation 
of a Pfaffian costs 0{M^) operations. 

B.2.2 Cayley's identity 

Cayley showed a useful identity:^^ 



det 



A12 .--Aim 
612 ... A2M 



= Pf 



X Pf 



A12 

-A12 



-612 . 
612 



■ Aim 

■ A2M 



A2M 



(B-IO) 

From this identity and the cofactor expansion of determi- 
nant, wc can obtain the relation between the Pfafhan and 
the inverse matrix of a skew-symmetric matrix A and the 
Pfaffian of a skew-symmetric matrix B which has same 
elements of A except for a-th row and column: 

detB detA^^A- 



Pf B 



1 -1 b 



Pf A 



Pf A 



m 

(B-11) 



where bm {m = 1, • • • . M) arc the updated elements of 
a-th row in B. eq. (B-11) allows us to calculate Pf B with 
0{M) operations from Pf A and A~^. 

B.3 Update technique for skew- symmetric matrix 

We derive one of the most important techniques for 
VMC with Pfaffians, which is similar to the update tech- 
nique for VMC with determinants. ^° 

B.3.1 Preparation 

For any nonsingular matrix A and any column vector 
u and V with the condition 1 + t)-'' Am 7^ 0, we have the 
Sherman-Morrison's formula:^^ 

[A + uv^] = A-^ - ^^^TA-i^ E Ar^UmVuA-l 

m,n 

(B-12) 

If we take A as a skew-symmetric matrix (A^ = —A), 
then we can derive an inverse matrix of B = A -|- uv^ — 

.T. 



vu 



(B-13) 

B.3. 2 Update formula for inverse matrix 

An inverse matrix of B with updated a-th row and 
column from the original skew-symmetric matrix A are 
calculated by 



4-1 



+ 



(B.14) 

where (m = 1, • • • , M) are the updated elements of 
a-th row in B and 5ij is the Kronecker's delta. If B is a 
singular matrix (Pf B = 0), the following formula holds 
instead of eq. (B-14): 



[PfB.B-i]^.=PfA 



4-1 



A-} + 



E^7™^" 



(B-15) 



where [Pf B • B ^] is a symbolic notation. Above formulas 
can be derived by using eq. (B-13) and by taking 



Ui 
Vi 



■ h- — A ■ 



(B.16) 



Equations (B-14) and (B-15) allows us to calculate B ^ 
and [Pf B - B'^] with 0(M2) operations from Pf A and 
A~^, which may be compared with 0{M^) operation if 
one calculates from scratch. 
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